Vessel delineation in magnetic resonance angiographic images

ABSTRACT

Delineating vessels in an angiogram involves two methods; graph generation and skeletonization. Generating a graph includes obtaining a digital image of an angiogram, recognizing a first growth point within the image, and identifying region boundary points around the growth point. The region boundary points are connected to the first growth point, thereby creating edges of a graph. The boundary point that has the greatest intensity is then selected as a second growth point, and additional region boundary points around the second growth point are identified. The additional region growth points are connected to the second growth point. The region boundary point with the greatest intensity in the image is then selected as a third growth point, and the method repeats until each point in the image is connected to another point in the graph. The skeletonization of the graph begins with recognizng a point in the graph as an endpoint of a vessel. This may be done explicitly through manual or automatic selection of specific points. It may also be done implicitly through a trimming process whereby graph branches of fewer than a certain number of connected points are discarded. The endpoints in the remaining branches are recognized as vessel endpoints. The skeletonization concludes with display of the delineated vessels. This may be done by superimposing the vessels in two or three dimensions over a conventional two-dimensional angiographic image such as a maximum intensity projection (MIP).

TECHNICAL FIELD

[0001] This invention relates generally to magnetic resonanceangiography (MRA). More particularly, the invention relates toprocessing magnetic resonance angiographic images, or angiograms, todelineate certain vessels in an angiogram.

BACKGROUND

[0002] Magnetic resonance angiography (MRA) is the magnetic resonanceimaging of the blood vessels in the body. In MRA, special pulsesequences are used by an MR scanner to cause flowing blood to appearvery bright and stationary tissue to appear very dark. If arterialstructures are being studied, additional pulses are applied to erase thesignal in veins. Multiple thin slices are obtained at adjacent levelsthrough the region of interest. In prior techniques, a computer thenstacks these images and creates a three-dimensional image. Theconstructed image can be rotated 360 degrees so that the vessels can bestudied in all projections.

[0003] MRA has become a primary method for the evaluation of vascularpathologies and increasingly for the purposes of surgical planning.Rapidly improving data acquisition methods have greatly improved imagequality. Much less attention, however, has been focused onpost-processing techniques for enhancing what is captured in the image.One such method, the maximum intensity projection (MIP), has become thestandard for vascular visualization. Because of its simplicity andnon-parametric basis as well as its high visual quality, the MIP isgenerally advocated when the angiogram is of high quality. However,limitations in image quality in MRA persist when it has been applied tomore challenging conditions such as in the abdomen, extremities, theheart, and where the vascular tree is highly overlapped such as in thecerebral MRA. The limited image quality affects the accuracy of analysisof vessel shape for diagnosis of vascular disease and the accuracy ofdetermination of vessel paths for surgical planning.

[0004] Specifically, there are several difficulties inherent in imageanalysis of the vasculature. One is the lack of definition in thevascular structure; within any finite resolution image, the distalextent of the vascular tree is indeterminate. Thus, the number andextent of detected vessels is dependent on acquisition and analysismethods. Another difficulty is that, even in moderately large vessels(by axial dimension), the vessel diameter will tend to be small relativeto the image resolution. A third difficulty is that vascular shape andanatomy is generally quite complex and variable.

[0005] An objective of the invention, therefore, is to provide a methodfor delineating vessels in angiograms. Another objective of theinvention is to provide such a method that works with vascular images ofthe type normally obtained from a magnetic resonance angiogram.

SUMMARY

[0006] Delineating vessels in an angiogram in accordance with theinvention involves two methods: graph generation and skeletonization.

[0007] Generating a graph includes obtaining a digital image of anangiogram, recognizing a first growth point within the image, andidentifying region boundary points around the growth point. The regionboundary points are connected to the first growth point, therebycreating a region of a graph. The boundary point that has the greatestintensity is then selected as a second growth point, and additionalregion boundary points around the second growth point are identified.The additional region growth points are connected to the second growthpoint. The region boundary point with the greatest intensity in theimage is then selected as a third growth point, and the method repeats.In one aspect of the invention, the method is repeated until each pointin the image is connected to another point in the graph.

[0008] The skeletonization of the graph begins with recognizing a pointin the graph as an endpoint of a vessel. This may be done explicitlythrough manual or automatic selection of specific points. It may also bedone implicitly through a trimming process whereby graph branches offewer than a certain number of connected points are discarded. Theendpoints in the remaining branches are recognized as vessel endpoints.

[0009] In one aspect of the invention, the method of graph generationmay be modified to handle circumstances where vessels overlap. In thiscase new region boundary points around the second and subsequent growthpoints are identified, wherein the new region boundary points areadditional region boundary points or present region boundary pointswhose distance within the graph from the present growth point is atleast a certain number of points. This introduces cycles into the graph.The skeletonization method is also modified, so that these cycles areidentified and then broken before the endpoints in the remainingbranches are recognized as vessel endpoints.

[0010] The skeletonization concludes with display of the delineatedvessels. This may be done by superimposing the vessels in two or threedimensions over a conventional two-dimensional angiographic image suchas a MIP. With the delineation provided by this method, a view of theimage can better detect the location and extent of vessels of interest.

[0011] Although the method has been developed initially for use inmagnetic resonance angiography, it is applicable to digital imagesproduced by other techniques, such as images generated with computedtomography angiography or x-ray angiography.

BRIEF DESCRIPTION OF THE DRAWINGS

[0012]FIG. 1 is a block diagram of a typical MRI scanner with which theinvention can be used.

[0013] FIGS. 2A-C illustrate a basic ORG method in accordance with theinvention for generating a graph representation of a low-resolutiondigital image and a method.

[0014]FIG. 3 shows a method for delineating vessels within an angiogramusing a graph constructed in accordance with the basic ORG method.

[0015] FIGS. 4A-B shows delineated vessels within several angiograms inaccordance with the delineating method.

[0016]FIG. 4 illustrates an alternative method for delineating vessels.

[0017]FIG. 5 shows the enhancement of an angiogram in accordance withthe alternative method.

[0018]FIG. 6 shows with several angiograms the effect of vessel overlapon the basic method.

[0019]FIG. 7 illustrates a problem that can occur in certainapplications of the basic method.

[0020]FIG. 8 shows the enhancement of an angiogram in accordance with amodified method that solves the problem noted in FIG. 7.

DETAILED DESCRIPTION

[0021]FIG. 1 is a block diagram of a typical MRI scanner 10. MRIscanners use the technique of nuclear magnetic resonance (NMR) to induceand detect a very weak radio frequency (RF) signal that is amanifestation of nuclear magnetism. The term “nuclear magnetism” refersto weak magnetic properties that are exhibited by some materials as aconsequence of the nuclear spin associated with their atomic nuclei. Inparticular the proton, which is the nucleus of the hydrogen atom,possesses a nonzero nuclear spin and is an excellent source of NMRsignals. The human body contains enormous numbers of hydrogenatoms—especially in water and lipid molecules. Although biologicallysignificant NMR signals can be obtained from other chemical elements inthe body, such as phosphorous and sodium, the great majority of NMRstudies are based on signals originating from protons that are presentin the lipid and water molecules within the patient's body.

[0022] The patient to be imaged is placed on a platform 12 and subjectedto several magnetic fields simultaneously or sequentially applied toelicit the desired NMR signal. The scanner 10 utilizes a strong staticmagnetic field in conjunction with gradient coils 14 and one or more RFcoils 16. The gradients and radio frequency components are switched onand off in a precisely timed pattern, or pulse sequence, by a pulsegenerator 18. Different pulse sequences are used to extract differenttypes of data from the patient. The NMR signal generated in response tothe coils is captured by a RF transceiver 20 and converted to digitalform by a data acquisition component 22. The scanner 10 is under thecontrol of a general-purpose computer 24. Operator commands are giventhrough an input device 26 such as a keyboard. The data acquired may bedisplayed as text and/or images on appropriate displays 28, 30,depending on the sophistication of the computer's software. Mass storage32 provides the computer 24 with software for generating and processingimages from acquired patient data and provides storage for that data.

[0023] The present invention is useful for processing magnetic resonanceangiographic images, or angiograms, generated by MRI scanners such asthe scanner 10. This processing delineates vessels otherwise not clearlyvisible to a viewer in an angiogram. In accordance with the invention,this processing can occur in real time, as the image is being generatedby the scanner, or it can occur later, using an electronic copy of theimage. The invention can be embodied in software that performsalgorithms in accordance with the invention. This software can be storedin mass storage 32 and loaded into the computer 24 or it can be storedand used on another computer that has access to an electronic copy ofthe angiogram.

[0024] It should be understood that the embodiment described herein isjust one example of how the invention can be used. It is not limited toenhancing MRA images. It may also be used for delineating vessels in anytype of image, including but not limited to images of the vasculatureobtained from computed tomography angiography or x-ray angiography.

BACKGROUND

[0025] Segmentation and Vessel Tracking Methods

[0026] There are two contrasting computational approaches to the problemof detection and characterization of vascular shape in MRA. The first isthat of segmentation or surface reconstruction which seeks to determinethe boundary between the interior and exterior of the vessels. Thisdetermination is inherently problematic because there are no definitivedistal endpoints of the vessels within the resolution of the images.Rather, the intensity of the vessels gradually diminishes as they becomesmaller due to the partial volume effect. Thus, intensity-based methodsof segmentation will always systematically underestimate the distalextent of the vessels and the size of the more distal vessels. Anotherproblem for direct segmentation, particularly of the smaller vessels, isthat points within the vessels are not interconnected to the extent thatthey are in more convex structures; relatively small imperfections inthe image may result in relatively large discontinuities in the vascularregion. Thus, region-based strategies for segmentation or surfacereconstruction of the vasculature will also tend to fail in the smallervessels.

[0027] An alternate approach to detection/characterization of vascularstructure in MRA focuses on the elongated nature of the vessels whichallows them to be represented by a line or series of lines which followthe centers of the vessels. These lines are useful as a framework foranatomical analysis of the vasculature and the identification ofanatomic variants. Moreover, they may provide the basis for measurementsof vessel diameter which is clinically relevant to the diagnosis ofatherosclerosis. In addition to these useful properties, the underlyingelongated, tubular structure of the vasculature has been found to beuseful a priori information for the detection of vessels. In severaldetection methods, the image intensity is assumed to be relativelyconstant along the central axis and to fall off sharply in twoorthogonal directions, at least for a given “scale” or degree ofblurring of the image which is quite plausible. Points which satisfythis criterion can be determined as well as paths along the centers ofvessels. Several related methods use 2^(nd) order partial derivatives todetermine vessel directionality in an analytic manner. Such methodsprovide both the path of the vessel and a measure of vessel diameter.However, the degree of user interaction can be quite high, presumablydue to problems at bifurcations where no single direction of the vesselcan be identified.

[0028] Skeletonization

[0029] Skeletonization methods address the question of the location ofthe central axes as well as the underlying topology of the object. Insuch skeletonization methods, which are typically applied tobinary-valued or pre-segmented images for purposes such as opticalcharacter recognition, the skeleton representation is intended toidentify key features of the object in the image as well as theinterconnectedness of those features. To assure such relevance in theskeletonization, topological constraints are typically an essentialfeature of the skeletonization processes. For example, in the skeletalrepresentation of a solid object no loops or holes should exist. Onesuch skeletonization method is the Skeleton by Influence Zones (SKIZ)which is applied to binary-valued images for reducing thick boundariesbetween various zones to a single-pixel width. In this case, thetopology of the object in the image is assumed to be purely web-like (intwo dimensions) with no discontinuous branches. Other skeletonizationsallow for branches with varying degrees of sensitivity. In addition,some methods have been developed for three-dimensional skeletonizationwhich encompass the more complex three-dimensional topologies. One suchmethod has been specifically applied to skeletonization of the MRA.

[0030] Such skeletonization methods are quite dependent on thesegmentation of the image and will suffer from any inadequacies in thesegmentation process. The methods of the invention seek to integrate, tosome extent, gray-scale detection methods with topology constraints forpurposes of both vessel detection and characterization. The underlyingprinciple for such a gray-scale skeletonization is ordered regiongrowing (ORG) connectivity.

[0031] Ordered Region Growing Connectivity

[0032] An ordered region growing (ORG) connectivity algorithm inaccordance with the invention permits gray-scale skeletonization ofmulti-dimensional images. Like Dijkstra's Shortest Path algorithm andminimum spanning tree algorithms, the ORG algorithm establishes adirected acyclic graph representation of the image. As such, unique“optimal” paths are established from a seed point to all other points inthe image based on gray-scale intensities. In the application ofDijkstra's algorithm to the determination of paths through images, forexample, the “cost” of a particular path might be either the cumulativeimage intensity or the cumulative inverse image intensity (depending onwhether, respectively, an intensity valley or ridge is preferred) andthe shortest path is that which has the lowest cost associated with it.The ORG algorithm, however, is not based on a cumulative cost and hasdistinct advantages over such path detection methods. One importantadvantage is that, except for the directionality of the ORG graph, thegraph is largely independent of the seed point location. A secondadvantage is that the ORG method has no bias towards spatially shorterpaths. Thus, for example, the path determined by a cumulative-costalgorithm along a curved vessel will be biased towards the side of thevessel at the inside of any given curve or the side with the smallerradius of curvature. In contrast, the path determined by the ORG methodwill be centered in the vessel. In terms of an overall tree structurerepresenting the vasculature, cumulative-cost methods tend to form thosetrees with relatively uniform branch lengths, suppressing the formationof lengthy branches. The ORG method. on the other hand, more trulyrepresents the length of each branch.

[0033] The ORG algorithm is defined as follows. For an N-dimensionalimage I: Z^(N)→R, let R_(n) represent the voxels that define the regionat iteration n. Let B_(n) represent the boundary of R_(n) at iteration nand B_(n) ⊂R_(n). R₀ may be any point or points in the image which willbe referred to as seed points. B₀ is the boundary of that region. G_(n)is then the set of growth points at the n^(th) iteration. Where:Neighbor refers to the set of immediate neighbors of a single point,either eight- or 26-neighbor for the three-dimensional case,(26-neighbor is used in this example), and Max refers to the singlepoint of maximum intensity of a set of points:

s _(n)=Max(B _(n))={xεB _(n) |∀yεB _(n) ,I(x)≧I(y)}  (1a)

G _(n)=Neighbor(s _(n))\R _(n)  (1b)

B _(n+1)=(B _(n) ∪G _(n))\s _(n)  (1c)

[0034] Certain advantages may be gained by the selection of meaningfulseed points such as those at the origin of a vascular tree of interestso as to establish the correct directionality of the skeleton althoughthe connectivity, described below, is largely independent of theseed-point placement with small-scale variability occurring mainly inthe vicinity of bifurcations.

[0035] The connectivity structure can be described in the terminology ofgraph theory. Connections, in that context. are the graph “edges.”Because each point is the offspring of a single other point, theconnectivity can be described as a map in terms of the set of all pairsof connected points, C.

E _(n+1) =E _(n) ∪{s _(n) −g _(n) |g _(n) εG _(n)}  (2)

[0036] And the directionality, which may be meaningful by itself, can beexpressed as a mapping from a given point to its parent P: Z^(N)→Z^(N)

P(g _(n))={s _(n) |g _(n) εG _(n)}  (3)

[0037] As note above, this connectivity structure, or graph, istree-like such that one and only one path exists between any two pointsin the image. This structure can also be characterized as acyclic orwithout loops.

[0038] This connectivity structure can be interpreted as a set of pathsbetween pairs of points in the region overgrown by the ORG region andthese paths themselves are optimal according to the “greatest minima”criteria with only small-scale discrepancies at bifurcations. Accordingto this criteria, the path within the ORG region has a point of minimumintensity which is greater than or equal to that of any alternative paththrough the image. In a sense then, the paths between any two points inthe image produced by the ORG algorithm are the most “continuous”relative to alternative paths. Under conditions where the intensity ofthe vessels tends to peak towards the center of the vessels, which isoften the case in MRA due to partial volume effects and reduced flowrates at the edges of the vessels, the ORG method will form pathsclosely aligned with the central axis of the vessels.

[0039] Implementation of the ORG algorithm is straightforward. Theprimary computational cost is due to determining the maximum point onthe boundary of the growing region, s_(n), which could, in principle,require re-scanning of a large portion of the image at each growthiteration. Thus, such a direct implementation could approach N²operations, for an image with N points. Such a computation would not befeasible for many images including 3D medical images where N>10⁷.However, provided a finite number of gray-scale intensities are allowed,finding the maximum can be carried out by a radix sorting process whichrequires on the order of N*log₂(M) operations (M is the number ofdiscrete gray-scale intensities). The portion of the ORG for storing andretrieving boundary points using a sorting mechanism is shown below. Forstoring a given boundary point, the intensity of the point is INTENSITY.The number of points stored below a given subdivsion L_(n,m) is given bythe variable OCCUPANCY. A pointer to the top of each of the M stacks isTOP and a pointer from each of the points in the stack to the next lowerpoint in the stack is LINK. A pointer to the point to be stored is NEWand a pointer to the point of maximum intensity in the M stacks (andtherefore in the ORG boundary) is MAXPOINT.

[0040] Set Intensity(L_(n,m))=m^(M−n) and set OCCUPANCY for all L_(n,m)to 0

[0041] Storage of a new boundary point, NEW:

[0042] Initialize n, m=0

[0043] While (m<M) {

[0044] OCCPANCY(L_(n,m))=OCCUPANCY(L_(n,m))+1

[0045] Is INTENSITY≧Intensity(L_(n+1, m+1))?

[0046] If so, m=m+1

[0047] if no, m=m

[0048] n=n+1}

[0049] LINK(NEW)=TOP(L_(n,m))

[0050] TOP(L_(n,m))=NEW

[0051] Retrieval of the point of maximum intensity on the boundary,MAXPOINT

[0052] Initialize n,m=0

[0053] If (OCCUPANCY(L_(0,0))>0) {

[0054] While (m<M) {

[0055] OCCUPANCY(L_(n,m))=OCCPANCY(L_(n,m))−1

[0056] Is OCCUPANCY(L_(n+1,m+1))>0?

[0057] If so, m=m+1

[0058] If no, m=m

[0059] n=n+1}

[0060] MAXPOINT=TOP(L_(n,m))

[0061] TOP(L_(n,m))=LINK(TOP(L_(n,m)))}

[0062] ORG skeletal lines from p₁ to p₂ with respect to a given origin,s_(o) of the ORG is denoted S_(ORG,so)(p₁, p₂) and is simply the set ofconnected points connected by the ORG graph from p₁ to p₂. Formally, theskeletal path S_(ORG,so)(p₁, p₂) is such that:

[0063] ∃(x₁, x₂, x₃, . . . , x_(n)) where x_(i)εS_(ORG,so)(p₁, p₂),

[0064] x₁=p₁, x_(n)=p₂ and

[0065] ∀iε(1, n), (x_(i−1), x_(i))εE_(final), (x_(i),x_(i+1))εE_(final),

[0066] FIGS. 2A-C illustrate the basic ORG method more conceptually. Theregion of the graph in FIG. 2A begins growing at a growth point, s_(n)(at origin of arrows) in the angriographic image and spreads outwards toneighboring pixels or G_(n) (known as region boundary points andindicated as light gray boxes). The growth point at each iteration(FIGS. 2B and C) is the point of maximum image intensity amongst allregion boundary points, B_(n) (dark gray boxes). The pattern of growthis recorded as a set of edges, E_(n), which is the accumulated set ofarrows from each iteration. The numbers in the boxes indicate the imageintensities. Points in the grown region not on the boundary areindicated by the black boxes.

[0067]FIG. 2D shows the ORG method applied to a low-resolutiontwo-dimensional digital angiographic image at the left, originating atthe point indicated by the arrow. The ORG acyclic graph is shown atright (the directionality of the graph is not shown). The pathsconnecting each point are optimal. A significant path is determined fromthis connectivity structure and the indication of additional points(arrows on right). In this way, centerlines of the vessels in the MRAcan be detected by the indication of a single seed point at the proximalend of a vascular tree and points at the distal ends of all terminalvessel segments.

[0068] In summary, generating the above graphic representation of anangiogram includes obtaining a digital image of the angiogram,recognizing a first growth point within the image, and identifyingregion boundary points around the growth point. The region boundarypoints are connected to the first growth point, thereby creating aregion of a graph. The boundary point that has the greatest intensity isthen selected as a second growth point, and additional region boundarypoints around the second growth point are identified. The additionalregion growth points are connected to the second growth point. Theregion boundary point with the greatest intensity in the image is thenselected as a third growth point, and the method repeats. The method isrepeated until each point in the image is connected to another point inthe graph.

[0069] Skeletonization by Explicit Selection

[0070] The Algorithm

[0071] The most straightforward application of the ORG method is theskeletonization of a given vessel path based on the specification of twoendpoints. This explicit selection can be done manually or doneautomatically. in accordance with given criteria. The graph at the rightof FIG. 2D illustrates the skeletonization method applied to an ORGgraph, with a seed point at the bottom and two selected endpoints at thetop. In general, for a given MRA, there will be one or more vessel treesof interest in which it is desirable to determine the paths and thebranching pattern. Construction of the vessel tree from vessel paths isrelatively straightforward because it requires the assembly of pathsfrom the seed point(s) to all given downstream or distal pointsidentified by the operator.

[0072] Conceptually and in implementation it is simplest to specify ameaningful seed point or points to initialize the ORG method. The seedpoint is generally chosen at the origins of the vessels of interest. Ifthe ORG path from a given distal point to a seed point is valid, it willhave a uniform upstream directionality. Algorithmically, to determinesuch a path, the chain of connectivities is simply followed in anupstream direction through connected points (which have associatedconnectivity information) until a seed point is encountered. Such pathscan be determined in real-time and thus the interaction can be done on atrial-and-error basis and failures of the algorithm or in pointselection can be undone and potentially corrected. Such failures arevisually obvious appearing as non-anatomic chaotic paths.

[0073] Symbolically, for all points i within an image I, a skeleton, S,can then be defined as all those points from which the set ofuser-defined peripheral points, PP, are descended as well as thoseperipheral points themselves with the connectivity of these pointsdetermined by the ORG algorithm. The set, Descendents (i), for a givenpoint i represents all points which can be traced upstream though thatgiven point, or, more formally, the descendents are all the points whichwould be disconnected from the rest of the connectivity graph if thatpoint and its associated connections were removed. Given that a Path (a,b) is the set of points within the image which are traversed within theORG connectivity graph going from a to b, including a and b, theDescendents could also be expressed as:

Descendents(i)={dεI|iεPath(d, P(i)), d≠i}  (4)

[0074] Then:

S={iεI|Descendents(i)∩PP≠Ø}∪PP  (5)

[0075] Determination of such a skeletal path, given the connectivity P,is simple and entails only the retracing of the links of P backwarduntil reaching the origin of the ORG path.

[0076] Moreover, it is possible to establish multiple branching treeswithin the images by simply supplying multiple seed points to the ORGalgorithm. Thus, there is a mechanism to compensate for cases, forexample, where regions of a vascular tree are disjoined or are onlyjoined through overly large vessels such as the aorta where the ORGconnectivity is poorly defined. For example, different seed points areneeded for skeletonizing the superior mesenteric artery and the hepaticartery although they are connected through the aorta.

[0077] An Application

[0078] Angiography is essential for preoperative planning prior tohepatic perfusion for tumor therapy and prior to liver transplantation.While digital subtraction angiography (DSA) has traditionally been usedfor such purposes, MRA is becoming an increasingly viable and lessinvasive approach. An example of the use of an abdominal MRA and thepotential improvement offered by skeletonization is discussed below.

[0079] Hepatic MRA is acquired with gadolinium enhancement fromintravenous bolus injection with the three-dimensional spoiled gradientecho, time of flight sequence. Skeletonization endpoints were choseninteractively within the liver from inspection of MIP images of limitedslice ranges where the presence and location of vessels was quite clear.In this case only the most obvious vessels were selected. Then, seedpoints were selected at central locations in the circulation fromamongst which only one will be associated with any given endpoint.

[0080]FIG. 3A shows on the left a maximum intensity projection (MIP)image of the anatomy of hepatic circulation (hepatic artery indicatedwith arrow, approximate boundary of liver highlighted with dotted line.The image on the right shows the skeletonization produced by explicitspecification of endpoints superimposed on the MIP showing the paths ofthe smaller vessels with greater clarity. In this example, all distalendpoints of the vessels are identified by the user while thecorresponding skeletal paths are determined and displayed in real-time.These endpoints may be identified either in individual slices or in MIPviews. Determination of the ORG connectivity of a 512×512×57 imagerequires several minutes on a Pentium II 500 MHz processor but oncecomputed, paths from each endpoint back to a seed point are determinedand superimposed on the MIP in real-time.

[0081]FIG. 3B shows an alternative display, a three-dimensional displayof the hepatic artery generated in accordance with the invention issuperimposed on the anatomy of hepatic circulation.

[0082] Skeletonization by Trimming

[0083] The Algorithm

[0084] A second method of skeletonization requires only a single seedpoint on a given vascular tree to be provided. In this method, theskeleton is formed by a trimming process whereby the central axes of thevessels are discriminated from trivial branches according to branchlength (defined in terms of number of connected points). The method isillustrated schematically in FIGS. 4A and B. FIG. 4A shows the ORG graphbefore trimming and FIG. 4B show the graph after trimming. At any branchpoint in the graph, all branches less than some minimum length arediscarded, except for the branch with the greatest downstream length.With this skeletonization method, all significant vessels are delineatedwithout the need to specify vessel endpoints. All that is requires isthe specification of minimum branch length.

[0085] This criterion can be defined fairly concisely with set notationand as an algorithm. Let R be the region overgrown by the ORG in whichall points in the region are greater than some given threshold value, asin a conventional intensity-based region growing. Let Siblings(i) be theset of all points with a common ORG parent:

Sibling(i)={jεR|P(j)=P(i)}  (6)

[0086] Let L be the connectivity distance from a given point i to itsfurthest descendent within the limited ORG, or, effectively, the branchlength:

[0087] Where the connectivity distance, D: Z^(N)×Z^(N)→Z, is just thenumber of connections between any two points within the connectivitygraph. (This definition can also apply to two given edges in the graph.)

[0088] L(i): Z^(N)→Z defined as

L(i)=Maximum({D(d, i)|dεDescendents(i)})  (7)

[0089] MIN be the minimum branch length which qualifies for abifurcation; one parameter of the skeletonization. Then S is the set ofpoints constituting the trimmed skeleton:

T={iεR|L(i)<MIN−1}∩{iεR|∃kεSiblings(i), L(i)<L(k)}  (8a)

D={Descendents(t)|tεT}∪T  (8b)

S=R\D  (8c)

[0090] An algorithm for producing such a skeletonization is given inTable 1: TABLE 1 1. Compute ORG connectivity of given region R from seedpoint 2. Initialize branch lengths L(i) of all points to 0 3. Untilreach steady state (no further changes in values of L for any point):for all i in R if L(i) > L(P(i)) − 1 then L(P(i)) = L(i) + 1 4.Initialize all points in a “skeleton” image, Skeleton, to 1 5. For all iin R: for all j siblings of point i if L(i) < L(j) and L(i) < Min thenSkeleton(i) = 0 6. Until reach a steady state (no further changes in anySkeleton) for all i in R if (Skeleton(Parent(i)) = 0) Skeleton(i) = 0

[0091] In addition to using the topology of the ORG for the purposes of“trimming” as just described, the topology can also be used to guide thegrowth of the ORG itself. Specifically, the ORG growth which until nowhas been assumed to be limited by some minimum threshold value, can belimited by a given number of significant bifurcations which has moreobjective meaning in MRA than any absolute image intensity value. Thiscan be done during the ORG graph growth by checking, at each growthiteration, if the growth has resulted in any new points qualifying asbifurcations.

[0092] The application of the method described above is limited mainlyto cases where pure threshold-based region growing would producereasonable results; large “leakage” in the region growing would likelycorrespond to false branches in the ORG skeletonization. Such conditionshave been found in some cerebral MRA's and the results of theapplication there will be discussed in detail below.

[0093]FIG. 5 shows how a region of the cerebral vasculature (middlecerebral artery tree) was skeletonized by the trimming process. The MIPof the original region image is shown at left and two-dimensionalprojection of the skeletonization is shown at right. Bifurcation pointsare indicated as the darker points in the projection. Direction of flowof the vessels, not shown in this visualization, is also obtained fromthe skeletonization process. This skeletonization requires only theidentification of a single seed point at a proximal location on thevascular tree, the nulling of the vessel in the upstream direction, andthe specification of the desired number of bifurcations and minimumbranch length. In this particular example. only one connectivity errorwas found in this result as indicated by the arrow.

[0094] Application to Circle of Willis MRA

[0095] Perhaps the site most suitable for a vascular skeletonizationmethod is in Circle of Willis (COW) MRA. COW images were acquired with astandard Spoiled Gradient Recalled Echo (SPGR), three-dimensionaltime-of-flight sequence. Skeletonization by trimming was applied to thegray-scale COW MRA in several arterial trees with discrimination ofvessels by the total number of bifurcations and by a minimum branchlength. Specifically, the skeletonizations were applied to basilarartery and branches and the two middle cerebral arteries. For eacharterial tree in each of three subjects, skeletonization was appliedwith both five and 10 bifurcations allowed with minimum branch lengthsof 15 and 25 voxel connectivity distance units (the length of the branchis considered as the number of voxels within it as an approximation).One seed point was specified at a proximal location on the vascular treebeing skeletonized along with the nulling of a rectangular regionimmediately upstream of the seed point to prevent undesiredskeletonization of the vessel in the upstream direction. The testing ofthe algorithm as reported in Table 2 was performed on the full imagesbut the speed of the algorithm is more accurately reflected in theapplication to a smaller cropped portion of the image wherenon-optimized tasks such as initialization are less significant;increase in the size of the image has essentially no effect on the speedof the core elements of this algorithm. In a 100×200×100 region of theCOW image which encompasses all visible vessels downstream of the MCA,the entire skeletonization requires less than 5 seconds on the MIPS10000 195 MHz processor.

[0096] Performance of the skeletonization was conducted on three MRA'swhich are normal in the region of interest which each included twomiddle cerebral artery (MCA) trees and one basilar artery (BA) tree.Skeletonization was performed on each vascular tree using the samecombinations of the number of bifurcations (five and 10) and minimumbranch length (15 and 25 voxel units) across the three MRA's primarilyto allow for comparability and pooling of the results amongst the threeimages. The capability of obtaining a “standardized” skeletonization ofa given cerebral MRA is itself of potential value, although the abilityto distinguish variability in image quality from variability in vascularstructure and function is unknown at this point. Small variation in theseed point placement generally only produces variation in theskeletonization in the immediate vicinity of the seed point because theunderlying ORG connectivity and directionality is relatively unchangedby small changes in the seed point location.

[0097] Thus, the skeletonization was rated simply according to thecorrectness of its connectivity. The detection or discrimination problemis not addressed by the trimming method and was not under considerationin this study; the vessels were effectively segmented by region growingwhich in this study was done conservatively such that nearly allovergrown regions would be part of the vasculature. Determination of thecorrectness of the skeletons was done by visual inspection of the MIPimage using limited slice ranges (“Subvolume MIPS”)and isosurfacevisualizations. In the evaluation of the skeletonization, each segmentof the vascular tree arising from a bifurcation (ten segments for caseof five bifurcations, 20 segments for case of ten bifurcations) wasrated as either correct or incorrect. Segments that did not belong totrue vessels or were incorrectly connected at their downstream originwere considered as errors.

[0098] The results of this study for the basic form of theskeletonization by trimming are shown in Table 2 and a sampleskeletonization is displayed in FIG. 5 for the middle cerebral arteryunder the processing specifications of 10 bifurcations and the branchlength of 15 voxel units. One error in all of the skeletonizations wasdue to a leakage of the skeleton into a non-vascular region, all othererrors were due to incorrect connectivities due to the incorrect joiningof vessels in close proximity to one another. Such errors, if occurringunder a lower specified number of bifurcation will necessarily occur athigher number of bifurcations. Thus, it is seen from the table thatalmost all such errors occur amongst the segments from the earlierbifurcations or amongst the largest vessels. A greater number of sucherrors occur in the arteries descended from the middle cerebral arterieswhere there is a greater apparent density of vessels than in thevascular tree originating from the basilar artery.

[0099] Overall, the skeletons produced by the trimming process werelargely accurate in their connectivity especially in the less densebasilar artery trees where less than one connectivity error, on averageoccurred. A higher rate of about three errors per tree did occur in theMCA arterial trees. These errors could almost entirely be attributed toa high degree of overlap of intensities of the larger vessels passingnearby to one another. However, because such errors occur in relativeisolation to one another, more sophisticated methods have largely beenable to correct these errors as will be discussed in the followingsections. TABLE 2 Study Number of Connectivity Conditions Samples Errors(average) MCA (5, 15) 6 2.3 ± 1.5 MCA (10, 15) 6 3.0 ± 1.7 MCA (5, 25) 62.1 ± 1.5 MCA (10, 25) 6 3.5 ± 1.4 basilar (5, 15) 3 0.0 ± 0.0 basilar(10, 15) 3 1.0 ± 1.7 basilar (5, 25) 3 0.0 ± 0.0 basilar (10, 25) 3 2.0± 1.7

[0100] Algorithm applied to middle cerebral arteries (MCA) and basilarartery in three MRA's. Numbers in parenthesis in “study condition”column refer to specified number of bifurcations and minimum branchlength (in voxel units) respectively. The samples were taken from threesubjects with two MCA and one basilar artery trees from each subject.

[0101] Modification of the ORG Algorithm

[0102] Generally, the constraint that connectivities amongst points inthe image be purely branching or acyclic in nature is suitable to theanalysis of the centerline of the vasculature; ideally those centerlinesare just a subset of the total connectivity tree which can be identifiedby explicit selection of endpoints on the vascular tree or by a trimmingprocess, as described above. However, this topological constraint may insome cases be too inflexible. This has been found to be the case wherethe density of vessels is very high and vessels appear to intersect toform apparent loops in the vasculature as tends to occur in the COW MRA(in addition to the real anatomical loops of the vasculature in thatregion). In these types of situations, the ORG method may establishincorrect connectivity between the distinct intersecting vessels, whichleads to the incorrect directionality of the graph beyond that point andwhich necessitates the disjoining of a true vessel segment because ofthe strictly acyclic nature of the ORG graph. This problem isillustrated in FIG. 6. Under such conditions it is reasonable to relaxthe acyclic constraint on the ORG topology so as to permit a limitednumber of cycles in the connectivity graph. A framework and specificmethods for doing so will be discussed below and in the followingsections. The resulting modified ORG algorithm, when applied to imagesof dense vasculature, will then be shown to be an adequate basis for theskeletonization.

[0103] A reasonable criteria for allowing a limited number of cyclesinto the connectivity graph is to limit the size of the cycles or, inother words, specify that any cycle in the connectivity graph must havea length greater than or equal to some given minimum size. Applied tothe ORG, the growth process would then proceed as before, but inaddition to growing into areas outside the growth region whichcorrespond to establishing acyclic connectivity, growth is alsopermitted onto previously grown points, provided that the size of thecycle formed is greater than some minimum size. Implementation of amethod with approximately this behavior is described in the followingsection.

[0104] At points where paths along vessel segments are disjoined in theORG graph, continuity along the vessel is usually obvious at the locallevel and clearly a greatest-minimum path and if only the local regionwere considered, the ORG would join points to either side of thediscontinuity. This is essentially the basis of the modified ORG method.While the modified method is not a technique for explicitly extendingthe graph, it does tend to “rejoin” such breaks within vessels in such away as is consistent with a local application of the method. This isqualitatively based on the observation that paths produced by an ORGmethod along a ridge in the intensity structure will be nearlyequivalent to the path of an ORG graph from two opposing origins on thatintensity ridge which includes their first collision point.

[0105] Characterization of the growth pattern associated with a “breakpoint” corresponds to that of two opposing ORG growths if only the localregion is considered, such as in FIG. 6. A two-dimensional test image onthe left simulates the effect of distinct vessels appearing to overlapone another due to their close proximity (arrowhead at left). The rightimage shows the region overgrown by the ORG method (originated at smallarrowhead) at one given point in the process. In this example, the ORGmethod incorrectly connects the vessels at the overlap point and breaksthe connectivity within a vessel (large arrowhead on right). However,these break points are usually relatively obvious and can be mended asdescribed hereafter. Note also that within a subregion of the image(dotted square on right) the ORG method produces the correctconnectivity along the vessel.

[0106] One problem with the above-modified ORG algorithm is that, oncecycles have been introduced into the graph, directionality within thegraph becomes ambiguous. As such, identification of paths between agiven pair of points requires traversing of the graph in a multitude ofdirections significantly complicating the computation.

[0107] However, while the graph is cyclic in nature, within a givenregion of the image and the corresponding graph a given pair of pointswill probably not be part of more than one cycle. Thus, the cycle sizedetermination can be made assuming there is only one path between agiven pair of points which can be determined by traversing the graph inan upstream direction.

[0108] The modification to the connectivity of the basic ORG method(from equation 2) is described below. Basically, a set of new edges, E′,is added to the set of edges E of the basic ORG. The qualifications forthe edges in E′ are described in (9a). Essentially, any cycle in the ORGgraph must have a size greater than some minimum value, CS_(min).Furthermore, to enforce the assumption that cycles in the graph occur inrelative isolation from one another, all the cyclic edges in E′ arerequired to be a distance of CS_(min) from one another.

[0109] Let a-s_(n) be an edge which connects the growth point s_(n) atany iteration to one of its neighbors, a, which is already in the grownregion R, (aεNeighbor(s_(n))\G_(n)) (G_(n) defined in (1b)). Let CS_(n):Z^(N)→Z be a function which determines, for any given edge added to theORG connectivity at a given iteration, what the size of theminimum-sized cycle it is associated with. “Arbitrary” is a functionthat arbitrarily chooses one element of a set that in this case isdetermined by the scan order at each iteration.

E _(pot,n) ={a−s _(n) |CS _(n)(a−s _(n))>CS _(min) }∩{a−s _(n) |∀e′εE′_(n) ,D(e′,a−s _(n))>CS _(min)}  (9a)

E′ _(n+1) =E′ _(n)∪{Arbitrary(E _(pot,n))|E _(pot,n)≠Ø}  (9b)

E _(n+1) =E _(n) ∪{g−s _(n) |gεG _(n) }∪E′ _(n+1)  (9c)

[0110] In summary, the method of graph generation may be modified tohandle circumstances where vessels overlap. In this case new regionboundary points around the second and subsequent growth points areidentified, wherein the new region boundary points are additional regionboundary points or present region boundary points whose distance withinthe graph from the present growth point is at least a certain number ofpoints. This introduces cycles into the graph, which are removed in themodified skeletonization method, as described below.

[0111] Modified Skeletonization by Trimming

[0112] The results of the basic skeletonization by trimming method,described above, can be improved using the modified ORG method such thatthe connectivity and directionality of the skeletonization is quitereliable. The method for extending the skeletonization by trimming isheuristic and is based on the observation that, within the higher-orderbranchings of the vasculature, bifurcations do not occur in closeproximity to one another and even in so-called trifurcations such as inthe popliteal artery, there is substantial separation between twodistinct bifurcations. This observation is significant in the context ofskeletonization because just such a formation tends to be produced whentwo distinct vessels nearly intersect and are incorrectly joined in theskeleton. This situation is illustrated in FIG. 7. The left image hastwo vessels that nearly intersect. Skeletonization based on the basicORG method may incorrectly connect the two vessels at that point. Ratherthan producing the correct connections of the vessel, atrifurcation-like structure may be formed in the skeleton consisting oftwo nearby bifurcations (two dark dots indicated by the arrow in theright image).

[0113] A modified trimming method seeks to remove or minimizetrifurcation-like formations in the skeletonization by taking advantageof the alternative connectivity or edges of the modified ORG graph. Morespecifically, a skeletonization by trimming is carried out as describedabove without regard for edges specific to the modified ORG graph. Then,essentially, a variety of alternative acyclic skeletonizations areformed by the inclusion of one or more of the cyclic edges produced bythe modified ORG and a corresponding removal of one of the originaledges. For each of these alternative skeletonizations, the one in whichall the bifurcations have the greatest degree of separation from oneanother is preferred.

[0114] Stepwise, the algorithm proceeds as follows: TABLE 3 1. Applymodified ORG algorithm. 2. Determine the basic skeletonization bytrimming. 3. Determine if the basic skeletonization by trimming can beimproved, in terms of reducing any “trifucations”, by the inclusion ofone or more of the “reconnection” edges from the modified ORG graph. Foreach “reconnection” edge: a. Add points to basic skeleton to completethe cycle associated with the reconnection edge. For purposes ofevaluating tri- furcations, trim off any trivial branches produced bythis step (less than 5 voxel units long). b. Determine the most likelytrifurcation point within the cycle (where the skeleton literallydivides in three or where there are two bifurcations in closesuccession). Then, remove the trifurcation: If the trifurcation iscomposed of two bifurcations which are both within the cycle, simplyremove the segment of the skeleton between the two bifurcations so as toobtain an acyclic skeleton. Otherwise, consider the trifurcation to bedetected-but-not-corrected. In that case, arbitrarily remove an edge inthe skeleton adjacent to the trifurcation which is within the cycle soas to obtain an acyclic skeleton.

[0115] Performance of the algorithm in Table 3 was evaluated andcompared with that of the algorithm in Table 1. For these tests, aslight difference in the implementation of the modified ORG method thatdescribed in (9) was used in which the distance between any two cyclicedges in E′ was considered to be the maximum of the distance from eitherof the two points to their mutual bifurcation point in the ORG graph asopposed to the sum of the distance from each of the two points to theirmutual bifurcation point. The cycle size, CS, was determined in asimilar way. This difference in implementation is unlikely to produceresults significantly different from those that would be produced by adirect implementation of the algorithm described by (9). For thesetests, the minimum cycle size, CS_(min) was set to be the same as thatof the minimum branch length.

[0116] The results were assessed only in the MCA artery region wherethere were a sizeable number connectivity and directionality errors intrimming of the basic ORG graph. Table 4 is a comparison of the modifiedSkeletonization by trimming algorithm with the basic skeletonization bytrimming algorithm (Table 2). Table 5 shows that the net result of themodified skeletonization by trimming was that less than one connectivityerror occurred within the vessels skeletonized under all conditions. Thecomplete computation of the skeleton, after the interactiveinitialization of the algorithm, took less than 10 seconds on100×200×100 region. TABLE 4 Correctly Remaining Added Added MissingStudy Number of Connections Connections Connections Conditions Samples(average) (average) (average) MCA (5, 15) 6 1.0 ± 0.6 0.5 ± 0.5 0.3 ±0.8 MCA (10, 15) 6 3.2 ± 1.6 1.0 ± 0.8 0.2 ± 0.4 MCA (5, 25) 6 1.5 ± 0.50.8 ± 0.7 0.2 ± 0.4 MCA (10, 25) 6 3.5 ± 1.0 1.2 ± 0.9 0.2 ± 0.4

[0117] TABLE 5 Undetected Detected but Connectivity Uncorrected StudyNumber of Errors Errors Conditions Samples (average) (average) MCA (5,15) 6 0.8 ± 1.2 0.2 ± 0.2 MCA (10, 15) 6 0.3 ± 0.8 0.2 ± 0.2 MCA (5, 25)6 0.7 ± 1.1 0.3 ± 0.3 MCA (10, 25) 6 0.7 ± 1.1 0.2 ± 0.2

[0118] Artery-Artery Separation

[0119] The high density and overlapping nature of the cerebral vasculartree makes visualization of individual vessels in MIP projectionsdifficult, particularly when vessels are highly enhanced by high flowrates even with careful selection of the MIP orientation and slicerange. This problem has been addressed by the application of shadedsurface display to provide depth information and by a method wherebyboth the surface and medial axes are detected allowing for selectivecolorization. In a likewise, but more powerful way, the cerebralskeletonization can be used to improve the quality of the visualization.The method suggested here is to effectively disentangle or separate thearterial sub-trees. Given that centerlines of vessels can be identifiedreliably by the methods described in the previous sections specificvascular sub-trees can easily be identified, based on, for example, thespecification of the origin a particular vessel. Points nearby to theskeleton downstream of the indicated point are set to zero intensity,Specifically, good results have been found zeroing out points within two26-neighbor dilation iterations of the downstream skeleton. The resultof the application to the a cerebral vascular tree is shown in FIG. 8.The portion of the vascular tree downstream of a given point on thevascular tree (arrow on left image) is nulled according to theconnectivity of the modified ORG skeletonization within two iterationsof a 27-neighbor dilations. Note a small section of an artery has beenunintentionally nulled (arrow in right image) due to the close proximityof distinct vessels. This method is useful for obtaining an unobstructedview of individual vessels in the cerebral MRA.

[0120] Once the desired skeleton of the vascular tree has been produced,this disentangling procedure can be conducted in a very interactivemanner. The non-optimized speed is less than five seconds on a100×200×100-cropped region.

[0121] The invention addresses the problems of visualization andcharacterization of small vessels in MRA which are inadequately resolvedby other methods. These problems include that of detection of the pathsof smaller vessels under conditions where regional or intensity-basedsegmentation methods are inadequate in the abdominal MRA and the problemof determination of the proper connectivity or anatomical relations ofthe vessels such as in the cerebral MRA. Like several other methods,these problems are addressed within the framework of gray-scaleoperations, but the novel principles of ORG connectivity applied to thisproblem incorporates global gray-scale intensity properties andtopological constraints into the formation of the centerlines both ofwhich are absent from previous methods. Specifically, paths are formedaccording to the greatest-minima property such that they follow thecenter of the vessels, to within the resolution of the image grid,provided that the image intensity of the vessels peak towards theircenter, which is typical for smaller vessels in MRA.

[0122] In implementation, the methods are practical and fast. All themethods are interactive in nature but small-scale variations in userinput have relatively minor effect on the resulting path determinations.Methods presented for the detection of vessels using ORG connectivity inthe abdominal MRA or in the smaller vessels of the cerebral MRA are themost intensely interactive requiring the specification of all distalendpoints. Improvement in the automation of this method within thecontext of the ORG is entirely possible. The paths may be incorporatedinto improved methods for determination of vessel diameter.

[0123] Having illustrated and described the principles of the inventionin an exemplary embodiment, it should be apparent to those skilled inthe art that the illustrative embodiment can be modified in arrangementand detail without departing from such principles. For example. many ofthe software aspects of the embodiment may be implemented in hardwareand many of the hardware aspects may be implemented in software. In viewof the many possible embodiments to which the principles of theinvention may be applied, it should be understood that the illustrativeembodiment is intended to teach these principles and is not intended tobe a limitation on the scope of the invention defined in the followingclaims. I therefore claim as our invention all that comes within thescope and spirit of these claims and their equivalents.

I claim:
 1. A computer-implemented method of generating a graphicrepresentation of an angiogram, the method comprising: a. obtaining adigital image of the angiogram; b. recognizing a first growth pointwithin the image; c. identifying region boundary points around the firstgrowth point; d. connecting the region boundary points to the firstgrowth point, thereby creating edges of a graph; e. selecting the regionboundary point with the greatest intensity as a second growth point; f.identifying additional region boundary points around the second growthpoint; g. connecting the additional region boundary points to the secondgrowth point, thereby growing the graph; h. selecting the regionboundary point with the greatest intensity in the image as a thirdgrowth point; and i. repeating steps f through h with the third andsubsequently identified growth points.
 2. The method of claim 1 whereinsteps f through h are repeated until each point in the digital image isconnected to another point in the graph.
 3. The method of claim 1including recognizing multiple first growth points.
 4. The method ofclaim 1 wherein the first growth point is a seed point provided by auser.
 5. The method of claim 1 wherein the digital image is threedimensional and identifying region boundary point comprises identifyinga set of adjacent points in three dimensions.
 6. The method of claim 1wherein the points in the image correspond to pixels or voxels.
 7. Themethod of claim 1 wherein the angiogram is obtained from magneticresonance angiography.
 8. The method of claim 1 wherein the angiogram isobtained from computed tomography angiography.
 9. The method of claim 1wherein the angiogram is obtained from x-ray angiography.
 10. A methodin accordance with claim 1 for delineating vessels within an angiogram,the method including: recognizing a point in the graph as an endpoint ofa vessel; and defining a path between the endpoint and the first growthpoint through connected points.
 11. The method of claim 10 wherein therecognizing step comprises recognizing a user-specified end point. 12.The method of claim 10 wherein the recognizing step comprises:discarding branches in the graph of a fewer than a certain number ofconnected points; and recognizing the last point in a remaining branchas an endpoint of a vessel.
 13. The method of claim 10 wherein thedefining step comprises superimposing the path on a maximum intensityprojection of the angiogram.
 14. The method of claim 10 wherein thedefining step comprises generating a three-dimensional display of thepath.
 15. A computer-readable medium on which is stored computerinstructions for performing the steps of claim
 1. 16. An apparatus forgenerating a graphic representation of an angiogram, comprising: meansfor obtaining a digital image of the angiogram; and means for:recognizing a first growth point within the image; identifying regionboundary points around the first growth point; connecting the regionboundary points to the first growth point, thereby creating edges of agraph; selecting the region boundary point with the greatest intensityas a second growth point; identifying additional region boundary pointsaround the second growth point; connecting the additional regionboundary points to the second growth point, thereby growing the graph;and selecting the region boundary point with the greatest intensity inthe image as a third growth point.
 17. A computer-implemented method ofgenerating a graphic representation of an angiogram, the methodcomprising: a. obtaining a digital image of the angiogram; b.recognizing a first growth point within the image; c. identifying regionboundary points around the first growth point; d. connecting the regionboundary points to the first growth point, thereby creating edges of agraph; e. selecting the region boundary point with the greatestintensity as a second growth point; f. identifying new region boundarypoints around the second growth point, wherein the new region boundarypoints are additional region boundary points or present region boundarypoints whose distance within the graph from the present growth point isat least a certain number of points; g. connecting the new regionboundary points to the second growth point. thereby growing the graph;h. selecting the region boundary point with the greatest intensity inthe image as a third growth point; and i. repeating steps f through hwith the third and subsequently identified growth points.
 18. A methodin accordance with claim 17 for delineating vessels within an angiogram,the method including: discarding branches in the graph of a fewer than acertain number of connected points; identifying cycles in the graph;breaking each cycle in the graph; recognizing the last point in aremaining branch as an endpoint of a vessel; and defining a path betweenthe endpoint and the first growth point through connected points. 19.The method of claim 18 wherein the defining step comprises superimposingthe path on a maximum intensity projection of the angiogram.
 20. Acomputer-readable medium on which is stored computer instructions forperforming the steps of claim 18.